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Abstract 

We discuss various methods used to investigate the geometric com- 
plexity of earthquakes and earthquake faults, based both on a point- 
source representation and the study of interrelations between earth- 
quake focal mechanisms. We briefly review the seismic moment tensor 
formalism and discuss in some detail the representation of double- 
couple (DC) earthquake sources by normalized quaternions. Non-DC 
earthquake sources like the CLVD focal mechanism are also consid- 
ered. We obtain the characterization of the earthquake complex source 
caused by summation of disoriented DC sources. We show that com- 
monly defined geometrical fault barriers correspond to the sources 
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without any CLVD component. We analyze the CMT global earth- 
quake catalog to examine whether the focal mechanism distribution 
suggests that the CLVD component is likely to be zero in tectonic 
earthquakes. Although some indications support this conjecture, we 
need more extensive and significantly more accurate data to answer 
this question fully. 

Keywords: Earthquake focal mechanism, double couple, CLVD, 
quaternion, geometric barriers, statistical analysis 
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1 Introduction 



Several statistical methods can be used to study the geometric complexity 
of an earthquake fault zone or fault system. Some are based on representing 
an earthquake as a point source, and the geometric complexity of a source 
reveals itself in a complex structure of a seismic moment tensor. Another 
method is to investigate the geometric complexity of the fault system as 
expressed in a set of earthquake locations and their focal mechanisms. 

1.1 Models for complex earthquake sources 
1.1.1 Phenomenological observations 

• 1. Geologic and geophysical studies of earthquake focal zones point to 
significant complexity in the rupture process (see, for example, King, 1983, 
1986). Moreover, almost any large earthquake is now analyzed in detail with 
its rupture history represented in time-space-focal mechanism maps which 
usually exhibit an intricate moment release. Although such results suggest 
that an earthquake focal zone is more complex than the standard planar 
model (Aki and Richards, 2002), phenomenological investigations cannot de- 
scribe the fault patterns appropriate for its quantitative modeling. Such a 
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description needs to be based on statistical treatment of the observed fault 
geometry. 

1.1.2 Point source solutions 

• 2. As a first approximation, the moment tensor is a second-rank matrix 
(Backus & Mulcahy, 1976a,b; Backus, 1977a,b; Kagan, 1987). The main com- 
ponent of the second-rank tensor is a double-couple (Burridge and Knopoff, 
1964). The presence of a significant non-double-couple or a CLVD component 
(Knopoff and Randall, 1970) is one measure of source complexity (SectionHJ). 
In this work we assume that for earthquakes the isotropic component equals 
zero. 

Julian et al. (1998) and Miller et al. (1998) discuss various physical mech- 
anisms thought to be responsible for the non-DC and, in particular, CLVD 
earthquake sources. They review many published papers that describe the 
registration of the CLVD component. Based on theoretical considerations, 
they indicate that most of the CLVD cases should be seen in geothermal 
and volcanic areas and quote observations confirming this. However, the 
complexity of an earthquake source can also produce a non-zero CLVD com- 
ponent in some tectonic earthquakes. Thus, the CLVD evaluation by modern 
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moment tensor solutions should, in principle, quantitatively characterize the 
complexity of the earthquake fault zone. Julian et al. (1998) and Miller et al. 
(1998) provide many published examples of a significant CLVD component 
present in tectonic earthquakes. 

However, recent work by Frohlich and Davis (1999), Kagan (2003), and 
Frohlich (2006, pp. 228-235) appears sceptical that present inversion tech- 
niques can obtain an accurate estimate of the CLVD component for tectonic 
earthquakes. Fig. 16 in Miller et al. (1998) also demonstrates the lack of cor- 
relation between the CLVD component values obtained from different earth- 
quake catalogs, indicating that these non-zero CLVD values may come from 
systematic effects. These results suggest that routinely determined CLVD- 
values would not reliably show the deviation of earthquake focal mechanisms 
from a standard DC model. 

• 3. Higher-rank point seismic moment tensors were introduced by Backus 
and Mulcahy (1976a,b) and Backus (1977a,b). Silver and Jordan (1983), 
Silver and Masuda (1985), Kagan (1987, 1988), McGuire et al. (2001), Chen 
et al. (2005) considered various aspects of higher-rank point seismic moment 
tensors. Kagan (1987) argued that the third-rank seismic moment tensor 
should show the complexity of the earthquake source, i.e., its difference from 
the standard planar rupture model. However, with the data available now, 
inversion results indicate that only the extent and directivity of the rupture 
can be obtained by analyzing higher-rank tensors. 
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1.1.3 Set of earthquake focal mechanisms 

• 4. Study of higher-rank correlation tensors (Kagan and Knopoff, 1985b; 
Kagan, 1992a). These correlation tensors more completely describe the inter- 
relation between focal mechanisms and their spatial distribution. However, 
their interpretation is still difficult and the low accuracy of available catalogs 
makes conclusions uncertain. Below we consider one of the invariants of the 
correlation tensor (tensor dot product for two arbitrary solutions). 

• 5. Investigation of the geometric relations between the double-couple earth- 
quake mechanism solutions (Kagan, 1982, 1990, 1991, 1992b, 2000, 2005). 
These papers studied the interrelation between the pairs of earthquake focal 
mechanisms by determining the 3-D angle between two DC solutions (Sec- 
tion and considered the angle distribution in space-time. It was shown 
that the angle is distributed according to the Cauchy law and increases with 
time and distance interval between earthquakes. However, as we will see 
later, the fault complexity depends not only on the pairwise distribution 
of 3-D rotation angles, but also on the distribution of rotation poles on a 
reference unit sphere. 
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1.2 Geometric complexity of earthquake faulting 

There is similarity between the geometric complexity of the earthquake fault 
zone and that of the earthquake fault system. Their comparability is the 
result of a general self-similarity of earthquake occurrence: earthquake rup- 
ture propagates over a complex fault pattern. This pattern is then seen in 
occurrence of aftershocks and other dependent events. Therefore, we may 
assume that the fault pattern complexity, when considered for small time 
intervals, would approach the geometric complexity of each earthquake rup- 
ture. Fractal features of the spatial distribution of earthquake hypocenters 
and epicenters (Kagan, 2007a) support this conjecture. 

An important condition which contributes to the complexity of the earth- 
quake fault system is the compatibility of elastic displacement: there should 
be no voids or material overlap created in fault configurations (Gabrielov et 
al, 1996). The authors developed a mathematical framework for calculating 
the kinematic and geometric incompatibility in tectonic block systems, both 
rigid and deformable. They concluded that due to geometric incompatibili- 
ties at fault junctions, new ruptures must be created to accommodate large 
plate tectonic deformations. 

Historically, earthquake fault mechanics was patterned along the lines of 
engineering fracture mechanics (see, for example, Anderson, 2005), where 
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the tensile fracture (formation of cracks) in materials and bodies is a ma- 
jor concern. Voids are created in tensile cracks; therefore the displacement 
incompatibility condition is not satisfied. However, earthquakes occur in 
Earth's interior, where a considerable lithostatic pressure should prevent the 
appearance of voids. 

Moreover, engineering concepts as opposed to physical theories generally 
cannot be transferred into a new scientific field without major experimental 
work. This may explain why fracture mechanics did not significantly improve 
our understanding of earthquake rupture process (Kagan, 2006). 

Similarly, another engineering science discipline, the study of friction, 
is based on experiments with man-made surfaces. Though widely used in 
earthquake mechanics, its contribution to the theoretical foundations of this 
science is still uncertain (ibid). 

Here we concentrate on analyzing the CLVD component of the seismic 
moment tensor. The CLVD is the simplest measure of earthquake focal 
zone complexity. Only in Section [5] do we briefly review the dot product for 
seismic moment tensors. This product is the first invariant of the second-rank 
correlation tensor (see Section [T. 1.3[) and can also be used to characterize the 
complexity of the focal mechanism orientation. Further detailed study needs 
to address the other two invariants of the correlation tensor. 

Many measurements of focal zone geometry indicate that a planar earth- 
quake fault is only a first approximation; rupture is usually non-planar. How- 
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ever, it is important to know whether the focal zone of a single earthquake or 
the fault systems of many earthquakes can be represented by a distribution of 
small dislocations with the DC mechanism. If no CLVD component is present 
in tectonic earthquakes, one degree of freedom for each rupture patch can be 
excluded with a great savings in representing earthquake rupture patterns. 

Frohlich et al. (1989) and Frohlich (1990) studied the CLVD distribution 
for earthquakes with the rupture slip along surfaces of revolution and found 
that a certain geometric slip pattern produces a significant CLVD compo- 
nent. However, such smooth surfaces of revolution are unlikely during a real 
earthquake rupture. Our results (Kagan, 1992, 2000, 2007a) rather suggest 
that both the fault rupture system and focal mechanisms, associated with 
the rupture, are non-smooth everywhere. They are controlled by fractal dis- 
tributions. 

Below, in Section [2] and Section [3] we review the seismic moment tensor 
formalism and quaternion representation of the DC focal mechanism. Sec- 
tion H] discusses the CLVD component of a focal mechanism and theoretically 
evaluates the component for some models of the earthquake composite source. 
We use the global CMT catalog to investigate statistically (Section [6]) the 
spatial pattern of earthquake focal zones and the distribution of focal mech- 
anism rotations to infer whether the CLVD component may appear as the 
result of fault complexity. In Section [7] we discuss the challenges in studying 
the focal mechanism pattern and the possibility of using local catalogs of 
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focal mechanisms for these investigations. 

2 Tensor invariants 



A seismic moment tensor can be represented as a symmetric 3x3 matrix 



m 



m n m 12 m 13 
m 2 i m 2 2 rn 3 2 
m 31 m 32 m 33 



(1) 



and as such it has six degrees of freedom. The moment tensor is considered 
to be traceless or deviatoric (Aki and Richards, 2002). Hence its number of 
degrees of freedom is reduced to five. 

The eigenvectors of matrix ([T]) are vectors 

t = [1,0,0]; 



[0, 1,0]; 



;o, o, r 



(2) 



For known eigenvectors t and p, the DC source tensor can be calculated 
as (Aki and Richards, 2002, Eq. 3.21) 

m ij — jj[niUj + rijUi] 

= fi[(ti+Pi) (tj ~Pj)+ (tj + Pj )(t i -p i )}/2, (3) 

where fi is a shear modulus, n is a normal to a fault plane, and u is a slip 
vector (see Fig. [T]). Therefore, if we know the orientation of two eigenvectors, 
the moment components can be calculated. 
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The invariants of the deviatoric seismic moment tensor tensor m can be 
calculated as (Kagan and Knopoff, 1985a) 

Ii = Tr [m] = mu + "7.22 + ^33 

= Ax + A 2 + A 3 = , (4) 

where 'Tr' is a trace of the tensor and A, are the eigenvalues of a moment 
tensor. The second invariant or the norm of the tensor 

h — — (mnm 2 2 + mnm 33 + m 22 m 33 ) + m\ 2 + m\ z + m\ z 

= -(AiA2 + AiA 3 + A 2 A 3 ). (5) 

For a traceless tensor (TJJ 

h = rn 2 l2 + m 2 13 + mj 3 + (m J u + m 2 22 + m 2 23 )/2 

= lEJ:< = (A? + A^ + A^)/2, (6) 
z i=i j=i 

(Jaeger and Cook, 1979, p. 33). The scalar seismic moment is 

M = s/T 2 . (7) 

To normalize the tensor we divide it by M 

m = m'/M. (8) 

In the rest of the paper, unless specifically indicated, we use only the nor- 
malized moment tensors. 
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The third invariant is a determinant of a tensor matrix 

I 3 = Det[m] = m n m 22 m 3 3 + 2m l2 m 13 m 2 3 

- (mnm| 3 + m 22 mj 3 + m 33 m\ 2 ) = AiA 2 A 3 . (9) 

For a double-couple (DC) earthquake source 

/ 3 = 0, (10) 

i.e., 

min[|Ai|, |A 2 |, |A 3 |] = 0. (11) 
Thus, the normalized DC moment tensor has 3 degrees of freedom. 

3 Quaternions and representation of double- 
couple earthquake mechanism and its rota- 
tion 

3.1 Quaternions 

3.1.1 Normalized quaternions 

Kagan (1982) proposed using normalized quaternions to represent earthquake 
double-couple focal mechanism orientation. Quaternions are widely used to 
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describe 3-D rotations [group SO (3)] in space satellite and airplane dynamics 
(Kuipers, 2002), in geodesy (Shen et al, 2006), and in simulations of virtual 
reality, robotics and automation (Kuffner, 2004; Yershova and LaValle, 2004; 
Hanson, 2005). There are matlab packages available to perform quaternion 
math calculation (MathWorks, 2006), as well as the quaternion description 
in matematica (MathWorld, 2004). 
The quaternion q is defined as 

q = <?o + 9i i + <?2j + <?3k . (12) 

The first quaternion's component (q ) is its scalar part, qi, q 2 , and q 3 are 
components of a 'pure' quaternion; the imaginary units i, j, k, obey the 
following multiplication rules 

i 2 = j 2 = k 2 = -1; 

i x j = -J x i = k; 
k x i = — i x k = j; 

jxk=-kxj = i. (13) 

From (fl3l) note that the multiplication of quaternions is not commutative, 
i.e., it depends on the order of multiplicands. Non-commutability is also a 
property of finite 3-D rotations. Thus, in general 

q"xq' ^ q'xq", (14) 
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i.e., we need to distinguish between the right- and left-multiplication. 
The conjugate q* and inverse q -1 of a quaternion are defined as 

q* = <?o - Qi'i - <?2j - <?3k and q x q 1 = 1. (15) 

The normalized quaternion q = [q , q±, q 2 , #3] contains four terms which 
can be interpreted as defining a 3-D sphere (S 3 ) in 4-D space: 

Qo + + Q2 + <7 3 2 = 1 • (16) 

Hence the total number of degrees of freedom for the normalized quaternion 
is 3. For the normalized quaternion 

q* = q- 1 • (17) 



3.1.2 Quaternions and 3-D rotations 

The normalized quaternion can be used to describe a 3-D rotation: in this 
case the first term in ffTBT) represents the angle of the rotation and the fol- 
lowing three terms characterize the direction of its axis (Kagan, 1991). 

We use normalized quaternions to calculate a rotated vector R(v) by 
applying the rules of quaternion multiplication (fl~3|) 

R{y) = q x v x q- 1 . (18) 
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The vector v = [91,92,93] is represented in (TT8T) as a pure quaternion, i.e., 
its scalar component is zero. In 0181) the quaternion q is a rotation operator 
and the pure quaternion v is an operand (Altmann, 1986, p. 16). Similarly 
to (fT8l) the whole coordinate system can be rotated (Kuipers, 2002). 

We use normalized quaternion multiplication to represent the 3-D rotation 
of the DC earthquake sources. The quaternion multiplication is 

s = q x r . (19) 

The above expression can be written in components (Klein, 1932, p. 61) 

so = 90^0 - qin - 92^2 - 93^3 ; 

si = 91^0 + 9o^i ± 92^3 =F 93^2 ; 

52 = 92^0 + 9o^2 ± 93^1 =F 9i r 2 ; 

53 = 93^0 + 9o^3 ± 9i r 2 =F 92^1 , (20) 

where the upper sign in ± and =F is taken for the right-multiplication and 
the lower sign for the left-multiplication: s = r x q. 

Kuipers (2002, p. 133) indicates that the right-multiplication corresponds 
to the 3-D rotation of an object, whereas the left-multiplication is the rotation 
of the coordinate system. Distinguishing these multiplications is especially 
important when considering a sequence of 3-D rotations 

q (n) = q' x q" X ... X q^ , (21) 
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is the right-multiplication sequence which we will use here. The correspond- 
ing rotation is anti-clockwise with the rotation pole located on a 2-D reference 
unit sphere (Altmann, 1986). 

3-D rotations for quaternions of opposite signs are equal 



where S is a transformation operator of a 3-D rotation corresponding to a 
quaternion q. This means that the group SO (3) of the 3-D rotations has a 
two-to-one relation to the normalized quaternions. Altmann (1986, Ch. 10) 
describes the complicated topology of rotations due to this representation. 

3.2 DC moment tensor and quaternions 

3.2.1 Quaternion representation of a single DC source 

Kagan (1982) represented the orientation of a DC source by a normalized 
quaternion. When applied to the DC parametrization, the identity quater- 
nion (zero rotation) 



S[q] 




(22) 



1 = [1, 0, 0,0], 



(23) 



is identified with the strike-slip DC source with plunge angles 



ctT = olp = 0°, and «b = 90° , 



(24) 
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and azimuths 



Pt = 0°, 



and p P = 90° , 



(25) 



(Kagan, 1991, 2005). Any other DC source corresponds to a quaternion 
describing the 3-D rotation from the reference DC source (Eqs. [231425]) . 

There are several possible representations of rotation in 3-D. Among 
the commonly used are Euler angles about coordinate axes (Kuipers, 2002, 
Ch. 4.3) and a rotation by the angle $ about a rotation axis. The rotation 
pole is the point where the rotation axis intersects a reference unit sphere. 
We use the latter convention in this paper since it is more convenient for 
the quaternion technique. For an arbitrary DC source, the value of the ro- 
tation angle and the spherical coordinates, 9 and 0, of the rotation pole on 
a reference 2-D sphere (S 2 ) are then 



where is an azimuth (360° > <p > 0°), measured clockwise from North; and 
9 is a colatitude (180° > 9 > 0°), 9 = 0° corresponds to the vector pointing 
down. 

We use the known correspondence between the orthogonal rotation matrix 



2 arccos (q ), 
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arccos [q^/ sin($/2)], 



arctan fa/qi), if < 0°, then <p = 360° + 4> 



(26) 



R 



h Pi h 
h p 3 h 



(27) 
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and the normalized quaternion (Moran, 1975, Eq. 6; Altmann, 1986, pp. 
52, 162; Kuipers, 2002, Eq. 5.11). We obtain the following formula for the 
rotation matrix 



R 



9o + 9? ~ 92 - 4 2 (-<?o<?3 + 9i 92) 2 (q q 2 + qiq 3 ) 
2 (9o93 + 9192) 9o - 9? + 9a - ?i 2 (-g 9i + 929s) 
2 (-9o92 + 9i9a) 2 (g 9i + 929s) 9o - 9? - 9a + 9a 



(28) 



to obtain the quaternion's components. The above formula can be obtained 
by applying f|T8|) to each of the original t, p, and b vectors ([2]). Kagan 
and Knopoff (1985, their Eq. 5) provide another expression for the rotation 
matrix using direction cosines of the axes [the first term in the second matrix 
row should be corrected as £m(l — cos$) — msin$]. 

Comparing fl27|) with fl28|) we derive quaternion components from the ro- 
tation matrix direction cosines (Kuipers, 2002, p. 169; Hanson, 2005, pp. 149- 
150). For example, if go is not close to zero 



g = (1/2)^+^ + 63 + 1; 

9i = (b 2 -p 3 )/(4g ) ; 

g 2 = (t 3 - &i)/(4g ) ; 

g 3 = (pi-t 2 )/(4g ), (29) 

Since as many as three of the quaternion components may be close to zero, it 
is computationally simpler to use the component with a maximum absolute 
value to calculate the three other components. 
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The DC seismic moment tensor in eigenvector coordinates is m = diag[ 1, —1, ]. 
For the general orientation of an earthquake focal mechanism, the seismic mo- 
ment tensor ([1]) can be calculated from the normalized quaternion as follows 
(Kagan and Jackson, 1994): 

mil = q\ ~ QqIqI ~ 2q\ Ql + 2 qIQq + 8<?i ^QsQo + 4, 

, o 2 2 r> 2 2 , 4 ^22, 4 

+ 2g 2 g 3 - 2q 2 q + q 3 - 6g 3 g + q ; 
"^12 = 4 (qfq 2 - qiql - q^q + q 3 qo) ] 
m 13 = 2 (qfq 3 - 3qlq 2 q - 3q x qlq 3 - q x q\ 

+ 3gig 3 go + iho + 3g 2 gf g - 92 ?o) ; 

m 2 2 = -q$ + Qqfqj - 2q\q\ + 2q\q\ + 9>qiq 2 q 3 q Q 

4,o22 o22 4,/-22 4 

m 23 = 2 (qlq + 3q\q 2 q 3 - 3gig 2 g + Sqiqlqo 

- qiql - qhs + q2ql - 3q 2 q 3 ql) ; 

m 33 = 4 (q\ql - q\ql - Aq^q^Q - q\q\ + q\ql) . (30) 
3.2.2 Mutual rotation of DC sources 

A more complicated algorithm is needed for the rotation of any DC source 
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into another. The methods of quaternion algebra can be used to evaluate the 
3-D rotation angle by which one DC source can be so transformed (Kagan, 
1991). Alternatively, the standard technique of orthogonal matrices can be 
applied to this calculation (Kagan, 2007b). 

Given the symmetry of the DC source (Kagan and Knopoff, 1985a; Ka- 
gan, 1990; 1991) the go term in ffT6l) can always be presented as the largest 
positive term in this parameterization. In particular, to obtain the stan- 
dard DC quaternion representation, we right-multiply an arbitrary normal- 
ized quaternion q by one of the elementary quaternions (Kagan, 1991): 

i= [0, 1, 0,0]; 
j = [0, 0, 1,0]; 

k= [0,0, 0,1], (31) 

if the second, third or fourth term has the largest absolute value, respectively. 
For example, for the largest second term, qx 

q" = q'xi. (32) 

If the resulting first term is negative, the sign of all terms should be reversed 
(see Eq. [22]). 

As the result of multiplication by expressions (13T1) the quaternion q be- 
comes 

q x i = [go, q±, <?2, 93]; 
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q X l = [q u -g , -<?3, <?2j; 



q x J = [92, gs, -go, 



q x k = [g 3 , -g 2 , gi, -go] ■ 



(33) 



The transformations (1221 [331 describe an eight-to-one correspondence be- 
tween an arbitrary normalized quaternion and a quaternion corresponding 
to a normalized seismic moment tensor. We call this operator Q. Therefore, 
the quaternion £ 



is a one-to-one quaternion representation of a DC focal mechanism. We call 
this operator Q. 

It is easy to check that all the eight quaternions of the Q operator (1341 
produce the same moment tensor (|30|) . Thus we can write 



where M ( ) is an operator (1301) converting a quaternion into a seismic moment 
tensor matrix. 

In principle, we can use the non-normalized variables. In this case the 
norm of a quaternion would correspond to that of the tensor (i.e., a scalar 
seismic moment M for a DC source, see Eq. [7]). However, the general de- 
viatoric tensor (Eq. HJ see also Section HJ) has five degrees of freedom even 
after it has been normalized. Hence it cannot be represented by a regular 
quaternion. 




(34) 




(35) 
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Similarly to (ITS]) rotation of a DC requires quaternion multiplication (as 
shown, for example in Eq. |2"T1) . The rotated DC source then needs to be 
converted into a DC standard quaternion representation using (I34I) . 

Thus, in our representation, an arbitrary quaternion is both a rotation 
operator and a DC source after simple transformations have been per- 
formed (Kagan, 2005). Although the quaternion does not have the advantage 
of clearly identifying the DC source properties, its benefits are obvious. Mul- 
tiple rotations of the DC source as well as the inverse problem determining 
the rotation from one source to another are easily computed using the meth- 
ods of quaternion algebra (Kagan, 1991; Ward, 1997; Kuipers, 2002). 

The Fortran program which determines the 3-D rotation of DC sources 
is available on the Web - 



ftp://minotaur.ess.ucla.edu/pub/kagan/dcrot.for (see also fortran90 adap 



tation of the programme by P. Bird http:/ /peterbird.name/oldFTP/2003107-esupp/Quaternion.f90 



Frohlich and Davis (1999) also discuss the program. Kagan (2007b) supplies 
simplified algorithms for calculating a 3-D DC rotation. These algorithms 
can be written in a few lines of computer code. 

3.3 Lambert azimuthal equal-area projection 

The DC focal mechanism has a symmetry of a rectangular box (Kagan, 
1991; 2005). Thus, it is convenient to use the Lambert azimuthal equal-area 
projection of an octant for a display of many distributions associated with 
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the DC source. The coordinates of the projection are (Snyder, 1987, p. 185) 
X = C cos(0') sin(0' -<#,); 

y = C [cos(^) sin(0') - sin(^) cos(0') cos(0' - 0' o ) ] ; 

G = / 2 

V 1 + sin(^) sin(0') + cos(^) cos(0') cos(0' - ' 

(36) 

where 180° > 0' > -180° is the centered longitude (c/. Eq. [26]), 90° > 0' > 
—90° is the latitude; 0' o and 0q are the coordinates of the projection center. 
For octant projection we use 4>' = 45° and 9' = arctan ^ ~ 35.26°. Then 
fl36|) can be simplified 

X = C cos(0') sin( </>"); 

y = -= [sin(0') - cos(0') cos(0") ] ; 
V3 



c 



^ (37) 



+ sin(6 |/ ) + \/2 cos(0') cos(0") ' 
where 0" = — 45°. Kagan (2005, Eqs. 26-29) provides an equivalent equal- 
area projection formula for an octant, if plunge angles ctj ff24l of a DC solution 
are known. 

4 CLVD and Gamma index 

The Gamma index (Kagan and Knopoff, 1985a; Frohlich, 1990; Richardson 
and Jordan, 2002) is 

r = T x ^' (38) 
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(see Eqs. I5HTU1). For a DC source 

T = 0, (39) 

(Eq. [TO]) . The T- index ranges from —1 to 1; |T| = 1 corresponds to a pure 
CLVD source (Knopoff and Randall, 1970; Kagan and Knopoff, 1985a). 

4.1 Gamma index for purely random rotation 

Kagan and Knopoff (1985a) considered the problem of a CLVD index distri- 
bution for a composite source 

N 

i=i 

where Rj is a random rotation matrix and Hj is its transpose. They showed 
that for the large number of summands (N) the resulting source has the 
uniform T-index distribution in the range 1.0 > T > —1.0. This result 
may be used to explain the non-zero T-index value, sometimes obtained for 
earthquakes with a complex fault zone, i.e., an earthquake source comprising 
several DC components of different orientation. However, there are both 
theoretical and observational arguments suggesting that the structure of a 
source is complex for tectonic earthquakes but precludes the appearance of 
a CLVD component. 

In a quaternion notation, (140]) can be expressed as 

N 

m ( q ) = M x r ) . ( 41 ) 

1=1 
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where the operator M is given by (j35p and r is a random rotation quaternion 
that can be obtained using Marsaglia's (1972) algorithm (see more in Kagan, 
2005). 

Kagan and Knopoff (1985a, their Fig. la) used simulation to show that 
for the sum of randomly oriented DC sources the T-index is distributed uni- 
formly over the interval [—1, 1] for a large number of summands. Even for 
two DCs, the distribution is close to uniform. This simplicity of the T-index 
distribution presents a significant benefit in characterizing the CLVD com- 
ponent. Many other measures of non-DC properties for an earthquake source 
have been proposed (see, for example, Julian et al. 1998, Eq. 18), but lack 
our statistical advantage. 

4.2 Gamma index for composite sources 

4.2.1 General considerations 

We will now explore the T-index distribution when the 3-D rotation is 
not completely random, i.e., if the rotation pole is preferentially located in a 
DC focal mechanism. For simplicity we assume that only two DCs comprise 
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the composite source (N = 2 in Eq. |4"T|) 

M(s) = q + xM(q x r) , (42) 

where for generality we assume that these DC sources have different scalar 
moments and their moment ratio is x = M'/ M" . Since the quaternion q 
is arbitrary, we can select it to be the identity (123]) . The general rotation 
quaternion is (see Eq. I2SJ) 

r = VT^A*; 

r\ = A sin(6 l ) cos(0) ; 

r 2 = A sin(^) sin(^) ; 

r 3 = A cos(^) , (43) 

where 

A = sin($/2), (44) 

(see Eq. EH). 

Using mathematica (Wolfram, 1999), we calculate the third invariant (Q 
for fl32J. We find that if the rotation axis is b, i.e., 6> = 0° in (1431). J 3 = for 
all $'s. Similarly, if the rotation axis is either the normal to the fault plane 
or is a slip vector 

4> = 45° or (f>= 135° and 9 = 90° , (45) 
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the invariant is also zero for all $'s. This means that the sum of two DC 
sources is again a DC for these rotations (Frohlich et ai, 1989). The same 
result (1 3 = 0) is obtained, if more than two DC sources are rotated around 
the b-axis or the n- and u-axes (see Eq. [3]) and then added with different 
weights (as in Eq. |4"2"1) . 

King (1986) specified two kinds of geometric barriers connected with a 
change of earthquake failure surfaces: conservative and non-conservative. 
The former structure does not require creating new faulting or void space. 
In our notation it would correspond to the second case (I45p . The first case is 
a non-conservative system: the incompatibility of displacement would require 
producing new ruptures (Gabrielov et ai, 1996). Therefore, with both kinds 
of barriers a complex geometric source is still a DC. Figs. [2] and [3] display 
cartoons of possible fault-plane and focal mechanism arrangements for these 
barriers. 

4.2.2 Small rotations 

If A in (jUJ) is small, we can keep only lower order terms in (j4"2]) and 
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obtain for the third invariant ([9]) of the sum 

J 3 (s) « 4 A 2 x (1 + x) sin 2 (9) f c os 2 (0) - sin 2 (0) 



x sin 2 (6 l ) cos(6 l ) sin(0) cos(0) . (46) 

The first term in the right-hand part of the equation suggests that the in- 
variant reaches maximum in the equatorial plane (8 = 90°) when the pole 
coincides with the t- or p-axes (0 = 0° or = 90°). Fig. H] displays the 
distribution of the T-index for the value of the rotation angle $ = 10° and 
for equal DC components (% = 1.0). 

In calculations we use an approximate formula (146]) ; the exact expression 
(obtained by mathematica) , which is too long to show here, yields almost the 
same answer. The T-index values in Fig. H] are very small; thus, a non-zero 
CLVD component is unlikely to be obtained in the moment solutions for an 
earthquake consisting of such subevents. 



4.2.3 Large rotations 



Formula fj46|) also suggests that if = 45°, for small values of the rotation 
angle $ and for non-equal DC components the third invariant is proportional 
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to sin 3 (<3>/2), i.e., it is close to zero. In particular, if x = 1, by employing = 
45° in (j46p . we obtain that 1% = for arbitrary values of the rotation angle $. 
This means that if the rotation pole is located on a fault- or auxiliary-plane, 
the sum of two focal mechanisms (original and rotated) has a zero CLVD 
component, i.e., it is a pure DC source. The effect can be seen in Fig. [51 
where the T-index, calculated with the exact formula, is displayed for equal 
DC components. In this case, if the rotation pole is close to the t- or p-axis, 
the resulting source is almost pure CLVD. 

This property can be demonstrated for a few simple arrangements of 
sources (see also Frohlich et al, 1989). For example, a sum of two DCs 
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can be rotated into s = diag[\/2, — a/2? 0]withJ 3 (s) = 0. This arrangement 
can be represented as two strike-slip events with a fault-plane rotated by 45°. 

Another example is the rotation by 120° around a pole with coordinates 
[l/v3, 1 / a/3, l/v3], a turn which exchanges the position of coordinate axes. 

s = diagfl, -1, 0] + diag[-l, 0, 1] = diag[0, -1, 1] . (48) 

Again I 3 (s) = 0. 

For rotation around t-axis, for example, 

s = diag[l, -1, 0] + diag[l, 0, -1] = diag[2, -1, -1] . (49) 
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Such a sum would have T = 1, i.e., it is a pure CLVD source (see Fig. E]). 

The results above demonstrate that if the rotation pole is located even 
randomly on nodal-planes, the resulting sum source is a DC. This may be 
relevant in searching for a non-DC component in seismic moment solutions. 
Frohlich and Davis (1999) and Kagan (2003) argued that the CLVD compo- 
nent is not presently measured with accuracy sufficient to study its properties. 

5 Tensor moment dot product 

The tensor dot product has been introduced by Kagan and Knopoff (1985b) 
as one of the correlation tensor invariants describing a complex earthquake 
fault pattern. Alberti (2006) used it to characterize the similarity of earth- 
quake focal mechanisms. We compare two methods for characterizing earth- 
quake source complexity: rotation angle $ and tensor dot product 



where the second source is rotated with regard to the first by the angle $. 

The dependence of the product D on 3 parameters of 3-D rotation f[2"61) 
calculated using mathematica, can be described as 



2 > D > -2, 



(50) 



D 




x -6 + QA 2 + A 2 cos 2 (6) + A 4 



sin 4 (9) 




(51) 
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where A is given by (|44p . 

In Figs. EH5]we display the D- value for three choices of the $ angle. The 
maximum change of the .D-value from D = 2, corresponding to identical DC 
sources, occurs when the rotation pole is at the b-axis. For $ = 90° rotation 
D = —2 at the b-axis and D = 1 for other axes (cf. Eq. EE9I) . Hence, the 
range of the D change is four times greater for the b-axis rotation than for 
rotation at other axes. Therefore, from the .D-value alone we cannot fully 
infer the complexity and coherence of DC sources. Higher order invariants of 
the correlation tensor need to be studied to better describe the fault pattern 
complexity (Kagan, 1992a). 

6 Focal mechanisms statistics 

The orientation of a DC source can be characterized by the following three 
quantities: a rotation angle ($) of the counterclockwise rotation from the 
first DC source to the second, and a location of a rotation pole on a refer- 
ence sphere f[2"6T) - colatitude, 9, and longitude, (ft (Kagan 1991, 2003). Thus, 
to fully study the distribution of earthquake focal mechanisms we need to 
investigate a six-dimensional manifold (fiber bundle): a product of the 3-D 
Euclidean space and the 3-D rotational distribution of a DC source. This 
presents a double difficulty - there is no effective procedure to display and 
study this pattern and we do not have sufficient data to evaluate the param- 
eters of the six-dimensional distribution. 
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Therefore, to better analyze the focal mechanism pattern, we should con- 
struct only marginal distributions of the full 6-D manifold. The following 
examples serve to illustrate the theoretical methods described in previous 
sections. 

In Section 16.21 we analyze the distribution of earthquake centroids pro- 
jected on a focal mechanism reference sphere for each earthquake. Section [6731 
presents the distribution of rotation angles between pairs of DC solutions, 
and Section 16.41 discusses the statistics of rotation poles. 

6.1 Earthquake catalog 

We study the earthquake distributions for the global catalog of moment ten- 
sor inversions compiled by the CMT group (Ekstrom et al, 2005, and refer- 
ences therein; see also references and earthquake statistics for 1977-1998 in 
Dziewonski et al, 1999). The catalog contains 26,865 solutions over a period 
from 1977/1/1 to 2007/3/31. Only shallow earthquakes (depth 0-70 km) are 
studied here. 

The CMT catalog includes seismic moment centroid times and locations 
as well as estimates of seismic moment tensor components (Dziewonski et al, 
1981; Dziewonski and Woodhouse, 1983). Each tensor is constrained to have 
zero trace (first invariant), i.e., no isotropic component. Double-couple (DC) 
solutions, i.e., with tensor determinant equal to zero, are supplied as well. 
Almost all earthquake parameters are accompanied by internal estimates of 
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error. 

From the original CMT catalog we created subcatalogs of well-constrained 
solutions. These datasets are obtained by removing from the catalog (a) the 
earthquakes lacking all 6 independent components of the moment tensor, (b) 
solutions with a large relative error, and (c) the solutions with large CLVD 
component (Frohlich and Davis, 1999; Kagan, 2000). Since the larger earth- 
quakes usually have smaller errors (Kagan, 2003), fewer of these events are 
removed. About 85% of m > 6.0 earthquakes are well-constrained, whereas 
for smaller earthquakes (m > 5.0) more than 2/3 have been deleted. 

6.2 Centroid distribution on a focal sphere 

Similarly to Table 1 in Kagan (1992b) in Table 1 here we show the distribu- 
tion of numbers of centroids in a coordinate system formed by the t-, p-, and 
b-axes of earthquake focal mechanisms. Each quadrant of the focal sphere 
is subdivided into 110 spherical triangles and quadrilaterals with equal area 
(consequently covering equal solid angles), and we calculate the number of 
times a centroid in projected into these cells. We normalize the numbers, so 
that the total number of the entries of Table 1 equals 11,000. Therefore, if 
the distribution of centroids in the tpb-system were completely random, all 
the numbers in the table would be equal to 100. The entries 'Col. Angle' give 
the colatitude angle corresponding to the lower edge of a segment consisting 
of one or several cells. The number of pairs for each segment is shown in 
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'Pair number' column. 

Table 1 shows a strong concentration of centroids near the plane going 
through the t— p and the b-vectors. Therefore this plane should correspond 
to the fault-plane and that going through the t+p and b-vectors should 
usually correspond to the auxiliary plane. Most of the strong earthquakes 
are concentrated in subduction zones. Fig. 3 in Kagan and Jackson (1994) 
illustrates that the fault-plane should include both the t— p and the b-vectors. 

Hue and Main (2003) suggest that the vertical errors in centroid deter- 
mination would strongly influence the direction statistics between the pairs 
of events. Even for M6 earthquakes, about 1/3 of shallow centroids are as- 
signed the depth of 15 km. This means that the depth could not be accurately 
evaluated. 

These results and those shown in Table 1 by Kagan (1992b) support the 
conventional model of an earthquake fault: a rupture propagating with slight 
deviations along a fault-plane. Figs. [3] and @] suggest that an earthquake 
occurring on such a fault system should have the CLVD component close 
to zero even if subevent focal mechanisms are significantly disoriented by 
rotation around the axis normal to the fault-plane (n-axis). 

6.3 Rotation angle statistics 

Fig. M displays cumulative distributions of the rotation angle $ for shal- 
low earthquake pairs separated by a distance of less than 50 km. We study 
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whether the rotation of focal mechanisms depends on where the second earth- 
quake of a pair is situated with regard to the first event. Thus we measure 
the rotation angle for centroids located in 30° cones around each principal 
axis (see curves, marked the t-, p-, and b-axes) of the first event. 

The curves in Fig. are narrowly clustered, and are obviously well approx- 
imated by the DC rotational Cauchy distribution (Kagan, 1990, 1992b). This 
distribution is characterized by a parameter k; a smaller K-value corresponds 
to the rotation angle $ concentrated closer to zero. Thus, all earthquakes, re- 
gardless of their spatial orientation, have focal mechanisms similar to a close 
event. Earthquakes in the cone around the b-axis correspond to a smaller 
K-value than the events near the other axes. These results are similar to the 
results shown in Fig. 6 by Kagan (1992b). 

In Fig. [10] we show how the distribution of the rotation angle $ depends 
on the magnitude threshold and the selection of the well-constrained earth- 
quakes in the catalog (Section 16.11) . It is clear from the diagram that 1) for 
larger earthquakes the angle distribution is generally more concentrated near 
zero, and 2) the well-constrained earthquakes also have a smaller angle be- 
tween pairs. These conclusions are easily explained by the higher accuracy of 
strong event solutions (Kagan, 2003): the disorientation of the DC solutions 
is caused mostly by solution errors. Kagan (2000) indicates that the rotation 
angle error is on the order 10° for the best solutions. Fig. H] suggests that if 
subevent constituents of an earthquake have small rotation (<3> < 10°), the 
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CLVD component would be close to zero. 

Below we show the pattern analysis for the most accurate solutions: the 
well-constrained earthquakes m > 6.0. The disadvantage here is their small 
number. 

6.4 Rotation pole statistics 

From Figs. 0HS1 one can infer that the rotation pole position strongly influ- 
ences the CLVD component of a complex source. In Fig. [TT] we show the 
distribution of the rotation poles for the second earthquake focal mechanism 
on a reference sphere of the first event. Because of the symmetry of the 
DC source, we reflect the point pattern at our reference sphere at the planes 
perpendicular to all axes. Thus, the distribution can be shown on an octant 
of a sphere. We use the Lambert azimuthal equal-area projection (1371) . 

For example, if in Fig. [11], a pole is shown near the b-axis, this would 
mean that in both mechanisms the axis has almost the same orientation: 
the second mechanism is rotated by an angle $ around an axis intersecting 
the sphere at the pole. The same pattern occurs for the poles that are close 
to other axes. In the diagram, the angle is 15° < $ < 30°. The points 
seem to concentrate between t- and p-axes. Fig. [12] shows the point density 
distribution, again confirming that the greatest pole concentration is located 
near a nodal-plane. 

In Table 2 we display the distributions of the rotation poles on a reference 
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sphere of the first event. In this Table, as in Table 1, we subdivide the positive 
octant of the sphere into 55 spherical triangles and quadrilaterals with equal 
area and then calculate the number of times the rotation axis intersects each 
of these cells. 

The upper chart in Table 2 shows the distribution of axes for the rotation 
angles of less than 15°. The distribution is randomly uniform as can be 
expected, because these rotations are caused, most probably, by random 
errors. For larger angles (the second chart, 15° < $ < 30°) the distribution 
shows that the stresses of the first earthquake influence the focal mechanism 
of the second earthquake in a pair. 

To illustrate the data and our technique in Table 3, we show the calcu- 
lations for 18 pairs of earthquakes in the lower-left corner cell of the second 
chart in Table 2. Earthquakes are distributed more or less uniformly at sub- 
duction zones. We can see again that many of them have been assigned the 
standard depth of 15 km. As explained above in this Subsection, the t-axes 
of both solutions are almost identical. Therefore, the rotation poles are con- 
centrated near the axis projection, both in Table 2 and in Fig. [TTJ Because 
the rotation angle is relatively small (15° < $ < 30°) the p-axes also have 
close orientation. 

However, the third and fourth charts of the Table calculated for the large 
angles of rotation display a different behaviour. The rotation poles are con- 
centrated near the b- and p-axes. Inspection of these earthquake pairs in- 
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dicates that 20 of 41 duos are concentrated in an area 200 by 200 km near 
the New Hebrides Islands - a region of complex subduction tectonics. These 
earthquakes likely result from the complicated slab geometry. 

We note that in Table 7 in Kagan (1992b) which is similar to Table 2, the 
plane orthogonal to the b-axis was incorrectly rotated by 90° and the results 
collapsed on an octant. Therefore, the number pattern at the former Table 
rows at p- and t-axis is a mixture of both distributions. 

7 Discussion 

We studied the distribution of the non-DC (CLVD) component for a com- 
posite earthquake source or a fault system. A source is considered complex 
if it consists of two or more events with a DC focal mechanism differently 
oriented. The theoretical computations detail conditions under which such 
a source would produce a non-zero CLVD component. For most of the geo- 
metric barriers proposed as common features in an earthquake fault system 
(King, 1983, 1986), the CLVD component should be zero. Frohlich et at. 
(1989) and Frohlich (1990) came to similar conclusions. 

We tried to statistically estimate whether the geometrical pattern leading 
to a zero CLVD could be confirmed by analyzing an earthquake catalog 
of focal mechanism solutions. However, the results shown in the previous 
Section are not fully convincing. They do, however, suggest limited support 
for this conjecture. 
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Why can't analysis of earthquake focal mechanisms in available catalogs 
definitively explain the CLVD component presence? To understand shear 
rupture properties, we need to investigate the behavior of the seismic moment 
tensor sums for distances between events close to zero. At such distances, 
the intrinsic geometrical conditions of the rupture would play a major role. 
At larger distances the slab geometry and its interaction with the upper 
mantle would significantly influence the focal mechanisms and pattern of 
earthquake hypocenters/centroids. Moreover, because of the low frequency 
waves used in the seismic moment inversion, centroids in the CMT catalog 
have low accuracy both in the horizontal plane and at depth. From Smith 
and Ekstrom's (1997) and Kagan's (2003) results, horizontal accuracy in 
a centroid location can be estimated on the order of 15-20 km. Centroid 
depth uncertainty should be the same or higher. In our computations we 
use distances between the events, which should increase the random error 
but decrease systematic uncertainties. The accuracy values quoted above 
are comparable with the slab thickness. This means that our distributions 
are influenced strongly by slab geometry. 

Brudzinski et al. (2007) studied slab geometry for the subduction zones 
using higher accuracy (less than 10 km, see above) hypocentral global earth- 
quake catalogs. They found that earthquakes concentrate at two layers: dou- 
ble Benioff zones or the upper and lower boundaries of a subducting slab. 
This means that large earthquakes which are listed in the CMT catalog are 
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mostly connected with the subducting plate geometry, not with any local 
geometry of the rupture zone. The geometry of deforming thin elastic sheets 
of material is a complex problem even if the sheet is deforming in free space 
(Marder et al, 2007). 

For earthquakes occurring at the double Benioff zones, the treatment of 
focal mechanisms should be modified because their symmetry properties dif- 
fer from those of earthquakes in the middle of a slab. For DC sources of 
earthquakes at the slab boundary, we should be able, in principle, to identify 
not only which of the focal planes is a fault-plane, but also to resolve the "up 
and down" of the fault-plane (Kagan, 1990, p. 576). Hence, the only sym- 
metry operation for such a focal mechanism is the identity fl23|) . In such a 
case, the quaternion visualization techniques described by Hanson (2005, Ch. 
21-23) can be implemented to picture both a slab surface geometry and the 
frames of earthquake rupture associated with the surface. However, the tech- 
niques proposed above are developed for smooth rotations. As we indicated 
earlier (Section 11.21) earthquakes occur on fractal sets and the orientation of 
their focal mechanisms is also controlled by a non-smooth fractal distribu- 
tion. Therefore, new methods need to be created for displaying earthquake 
focal mechanisms and their spatial pattern. 

While investigating earthquake spatial distribution (Kagan, 2007a), we 
were able in some degree to make a smooth transition from a 3-D distri- 
bution of hypocenters for small distances to a 2-D epicenter distribution for 
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distances comparable or larger than the thickness of seismogenic crust. After 
appropriate adjustments, the obtained values of fractal dimension are com- 
parable for both cases. It seems, however, that the methods employed for 
simple point distributions in 2-D and 3-D are much more difficult to imple- 
ment for tensor quantities. The challenges are two-fold: theoretical problems 
of displaying and interpreting the 6-D distributions and the more practical 
problem of no high quality, massive datasets of focal mechanism solutions. 

Several extensive local catalogs of focal mechanism solutions have been 
compiled (Pasyanos et al, 1996; Kubo et al, 2002; Hardebeck, 2006; Clinton 
et al, 2006; Matsumoto et al, 2006). They contain thousands of events. 
Many of these occur near local seismographic stations and, therefore, their 
solutions have good depth control. However, the largest catalogs are based on 
first-motion analysis and their solutions have a significantly lower accuracy 
than do catalogs of the moment-tensor inversions (Kagan, 2002, 2003). The 
latter catalogs are not yet sufficiently extensive, well-documented, or tested 
to be used for rigorous statistical studies. However, it is possible that in a 
few years these catalogs will be improved. 
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Figure 1: Schematic diagrams of earthquake focal mechanism. Equal- area 
projection (Aki and Richards, 2002, p. 110) of quadrupole radiation patterns. 
The null (b) axis is orthogonal to the t- and p-axes, or it is located on the 
intersection of fault and auxiliary planes, i.e., perpendicular to the paper 
sheet in this display. 
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Figure 2: Schematic diagrams of earthquake focal zone. Lines show a complex 
fault plane; focal mechanisms of earthquake(s) are displayed for each fault 
segments. 

(a) Rotation around b-axis, non-conservative geometric barriers. 

(b) Rotation around u-axis (slip vector), conservative geometric barriers. 
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(c) 
(d) 
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Figure 3: Schematic diagrams of earthquake focal zone. Rotation around 
n-axis (fault-plane normal vector), conservative geometric barrier. 

(c) View from side. 

(d) View from above. 
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Rotation and summation of DC sources, r distribution 
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Figure 4: T-index distribution for small rotation angles ($ = 10°). See 
Eq. HQ Two equal DC sources, one rotated compared to another. The axis 
angles are shown at octant equal-area projection fl37j) ; see also (Kagan, 2005). 
Dashed lines are boundaries between different focal mechanisms. Plunge 
angles 30° and 60° for all mechanisms are shown by thin solid lines. 




Figure 5: T-index distribution for large rotation angles (<& = 90°). Two 
equal DC sources, one rotated compared to another. Octant projection and 
auxiliary lines are the same as in Fig. HI 
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Figure 6: Distribution of the tensor dot product for rotation angle $ = 10°. 
Two equal DC sources, one rotated compared to another. Octant projection 
and auxiliary lines are the same as in Fig. HI 
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Figure 7: Distribution of the tensor dot product for rotation angle $ = 45°. 
Two equal DC sources, one rotated compared to another. Octant projection 
and auxiliary lines are the same as in Fig. HI 
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Rotation of DC source, tensor dot product 
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Figure 8: Distribution of the tensor dot product for rotation angle $ = 90°. 
Two equal DC sources, one rotated compared to another. Octant projection 
and auxiliary lines are the same as in Fig. HI The D- value is —2 at the b-axis 
and +1 at the t- and p-axes. 



59 



CMT 77-2007/03/31 , shal-R-0— 50km, M5.0: — all, : T, -. P, * B, k- Cauchy (k=0.1), m- random 
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Figure 9: Distributions of rotation angles for pairs of focal mechanisms 
of shallow earthquakes (depth 0-70 km) in the CMT catalog 1977/01/01- 
2007/03/31; centroids are separated by distances between 0-50 km, magni- 
tude threshold m w = 5.0. Circles - all centroids ; crosses - centroids in 30° 
cones around the t-axis; plusses - centroids in 30° cones around the p-axis; 
stars - centroids in 30° cones around the b-axis. Left solid line is for the 
Cauchy rotation with k = 0.1; right solid line is for the random rotation. 
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Figure 10: Distributions of rotation angles for pairs of focal mechanisms 
of shallow earthquakes in the CMT catalog; centroids are separated by 
distances between 0-50 km; centroids in 30° cones around the b-axis. 
Stars - CMT well-constrained earthquakes, magnitude threshold m t = 6.0, 
1977/01/01-2005/01/01; plusses - CMT catalog, magnitude threshold m t = 
6.0, 1977/01/01-2007/03/31; circles - CMT catalog, magnitude thresh- 
old m t = 5.6, 1977/01/01-2007/03/31; crosses - CMT catalog, magni- 
tude threshold m t = 5.0, 1977/01/01-2007/03/31; diamonds - CMT not 
well-constrained earthquakes, magnitude threshold m t = 5.0, 1977/01/01- 
2005/01/01. Left dashed line is for the Cauchy rotation with k = 0.2; right 
solid line is for the random rotation. 
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Figure 11: Distributions of rotation poles for pairs of focal mechanisms of 
shallow well-constrained earthquakes in the CMT 1977/01/01-2005/01/01 
catalog. Centroids are separated by distances between 0-50 km; magnitude 
threshold mt = 5.6; the rotation angle 15° < $ < 30°. 
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Figure 12: Map of rotation poles for pairs of focal mechanisms of shallow 
well-constrained earthquakes in the CMT catalog, see Fig. [TTJ 
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